ATAC and RNA exploration

here::set_here()

## File .here already exists in /Users/koenvandenberge/Dropbox/research/brainStat/hbcRegenGithub/scripts

suppressPackageStartupMessages({
  library(slingshot)
  library(SingleCellExperiment)
  library(cowplot)
  library(rgl)
  library(clusterExperiment)
  library(RColorBrewer)
  library(ggplot2)
  library(pheatmap)
  library(wesanderson)
  library(gridExtra)
  library(irlba)
  library(tidyverse)
})

Import data

sds <- readRDS("../data/finalTrajectory/sling.rds")
counts <- readRDS("../data/finalTrajectory/counts_noResp_noMV.rds")
# counts <- round(counts)
sce <- readRDS("../data/finalTrajectory/sce_tradeSeq20200904.rds")
load("../data/ALL_TF.Rda")
clDatta <- readRDS("../data/finalTrajectory/dattaCl_noResp_noMV.rds")
timePoint <- readRDS("../data/finalTrajectory/timePoint_noResp_noMV.rds")

Select activated and renewed HBC clusters

I will select

• activated HBCs as cells sampled in the 24h timepoint that are in the starting cluster
• renewed (‘resting’) HBCs as cells sampled in the 14D timepoint that are in the HBC cluster

cols <- c("#E41A1C", "#377EB8", "#4DAF4A", "#984EA3", "#FF7F00", "#FFFF33", "#A65628", "#F781BF", "#999999")
view3d( theta = 10, phi = 10)
# Datta labels
plot3d(reducedDim(sds), aspect = 'iso', col=cols[clDatta], alpha=.6, xlab="UMAP1", ylab="UMAP2", zlab="UMAP3", box=FALSE, top=TRUE)
plot3d(sds, add=TRUE, lwd=3, col="black")

# timepoints
plot3d(reducedDim(sds), aspect = 'iso', col=cols[timePoint], alpha=.6, xlab="UMAP1", ylab="UMAP2", zlab="UMAP3", box=FALSE, top=TRUE)
plot3d(sds, add=TRUE, lwd=3, col="black")

# cluster labels used for TI
clSds <- apply(slingClusterLabels(sds), 1, which.max)
plot3d(reducedDim(sds), aspect = 'iso', col=cols[clSds], alpha=.6, xlab="UMAP1", ylab="UMAP2", zlab="UMAP3", box=FALSE, top=TRUE)
plot3d(sds, add=TRUE, lwd=3, col="black")


actHBCId <- which(clSds == 8 & timePoint == "24HPI")
restHBCId <- which(clDatta == "HBC" & timePoint == "14DPI")
grpHBC <- rep(NA, length(c(actHBCId, restHBCId)))
grpHBC[1:length(actHBCId)] <- "activated"
grpHBC[is.na(grpHBC)] <- "resting"
grpHBC <- factor(grpHBC)
countHBC <- counts[, c(actHBCId, restHBCId)]

Dimensionality reduction

library(scran)
library(uwot)
library(scater)

## 
## Attaching package: 'scater'

## The following object is masked from 'package:clusterExperiment':
## 
##     plotHeatmap

logFracs <- log(sweep(countHBC+1e-5, 2, colSums(countHBC), "/"))
gv <- modelGeneVar(countHBC)
hvg <- getTopHVGs(gv, n=1000)
pc10 <- irlba::irlba(logFracs[hvg,], nv=10)
um10 <- uwot::umap(pc10$v %*% diag(pc10$d))

plot(um10, col=alpha(c("orange", "darkseagreen3")[as.numeric(grpHBC)], .75), pch=16, cex=1/2, bty='l')
legend("topright", c("activated", "resting"), pch=16, col=c("orange", "darkseagreen3"), bty="n")

Batch correction with Harmony

library(harmony)

## Loading required package: Rcpp

harmony_embeddings <- HarmonyMatrix(logFracs, meta_data=grpHBC, do_pca=TRUE)

## Harmony 1/10

## Harmony 2/10

## Harmony converged after 2 iterations

umHarmony <- uwot::umap(harmony_embeddings)

plot(umHarmony, col=alpha(c("orange", "darkseagreen3")[as.numeric(grpHBC)], .75), pch=16, cex=1/2, bty='l')
legend("topright", c("activated", "resting"), pch=16, col=c("orange", "darkseagreen3"), bty="n")

Integration: Processing and dimensionality reduction of each dataset separately

library(Signac)
library(Seurat)

## 
## Attaching package: 'Seurat'

## The following object is masked from 'package:SummarizedExperiment':
## 
##     Assays

library(patchwork)

## 
## Attaching package: 'patchwork'

## The following object is masked from 'package:cowplot':
## 
##     align_plots

library(ArchR)

## Loading required package: data.table

## 
## Attaching package: 'data.table'

## The following objects are masked from 'package:dplyr':
## 
##     between, first, last

## The following object is masked from 'package:purrr':
## 
##     transpose

## The following object is masked from 'package:SummarizedExperiment':
## 
##     shift

## The following object is masked from 'package:GenomicRanges':
## 
##     shift

## The following object is masked from 'package:IRanges':
## 
##     shift

## The following objects are masked from 'package:S4Vectors':
## 
##     first, second

## Loading required package: rhdf5

## Loading required package: magrittr

## 
## Attaching package: 'magrittr'

## The following object is masked from 'package:purrr':
## 
##     set_names

## The following object is masked from 'package:tidyr':
## 
##     extract

set.seed(1234)
addArchRGenome("mm10")

## Setting default genome to Mm10.

# scATAC-seq importing and preprocessing
# from 20200610_scATAC_injury_archr_samples_pairwiseHBC_v3.Rmd
archProj <- loadArchRProject("../data/scATAC/archR_samples_pairwise_v2")

## Successfully loaded ArchRProject!

## 
##                                                    / |
##                                                  /    \
##             .                                  /      |.
##             \\\                              /        |.
##               \\\                          /           `|.
##                 \\\                      /              |.
##                   \                    /                |\
##                   \\#####\           /                  ||
##                 ==###########>      /                   ||
##                  \\##==......\    /                     ||
##             ______ =       =|__ /__                     ||      \\\
##         ,--' ,----`-,__ ___/'  --,-`-===================##========>
##        \               '        ##_______ _____ ,--,__,=##,__   ///
##         ,    __==    ___,-,__,--'#'  ==='      `-'    | ##,-/
##         -,____,---'       \\####\\________________,--\\_##,/
##            ___      .______        ______  __    __  .______      
##           /   \     |   _  \      /      ||  |  |  | |   _  \     
##          /  ^  \    |  |_)  |    |  ,----'|  |__|  | |  |_)  |    
##         /  /_\  \   |      /     |  |     |   __   | |      /     
##        /  _____  \  |  |\  \\___ |  `----.|  |  |  | |  |\  \\___.
##       /__/     \__\ | _| `._____| \______||__|  |__| | _| `._____|
## 

clHBCMerged <- archProj$clHBCMerged
plotEmbedding(ArchRProj = archProj, colorBy = "cellColData", 
              name = "clHBCMerged", embedding = "UMAP_cor")

## ArchR logging to : ArchRLogs/ArchR-plotEmbedding-49602e03504b-Date-2021-11-04_Time-15-22-57.log
## If there is an issue, please report to github with logFile!

## Getting UMAP Embedding

## ColorBy = cellColData

## Plotting Embedding

## 1 
## ArchR logging successful to : ArchRLogs/ArchR-plotEmbedding-49602e03504b-Date-2021-11-04_Time-15-22-57.log

archrUmap <- ArchR::getEmbedding(archProj, embedding = "UMAP_cor")
colnames(archrUmap) <- c("UMAP1", "UMAP2")
plotEmbedding(ArchRProj = archProj, colorBy = "cellColData", 
              name = "treat", embedding = "UMAP_cor")

## ArchR logging to : ArchRLogs/ArchR-plotEmbedding-49605b8c8386-Date-2021-11-04_Time-15-23-00.log
## If there is an issue, please report to github with logFile!
## Getting UMAP Embedding
## ColorBy = cellColData
## Plotting Embedding
## 1 
## ArchR logging successful to : ArchRLogs/ArchR-plotEmbedding-49605b8c8386-Date-2021-11-04_Time-15-23-00.log

clHBCMerged2 <- as.character(getCellColData(archProj, "clHBCMerged")$clHBCMerged)
clHBCMerged2[clHBCMerged2 == "C1_Inj"] <- "Activated"
clHBCMerged2[clHBCMerged2 == "C2_Hybrid"] <- "Hybrid"
clHBCMerged2[clHBCMerged2 == "C3_UI"] <- "Resting"
archrUmap$clHBCMerged2 <- factor(clHBCMerged2)
archrUmap$treatment <- getCellColData(archProj, "treat")$treat
peaks <- archProj@peakSet
names(peaks) <- paste0("peak",1:length(peaks))
atacPeakMat <- readRDS("../data/scATAC/peakMatArchrHbc_v2.rds")
rownames(atacPeakMat) <- names(peaks)
geneActivity <- readRDS("../data/scATAC/geneMatArchrHbc.rds")
atac <- Seurat::CreateSeuratObject(counts = atacPeakMat,
                           assay = "ATAC",
                           project = "10x_ATAC")
atac[["ACTIVITY"]] <- CreateAssayObject(counts = geneActivity)
atac$tech <- "atac"
atac$treatment <- archProj@cellColData$treat
atac$clHBCMerged <- clHBCMerged
# process gene activity
DefaultAssay(atac) <- "ACTIVITY"
atac <- FindVariableFeatures(atac)
atac <- NormalizeData(atac)
atac <- ScaleData(atac)

## Centering and scaling data matrix

# process peak matrix
DefaultAssay(atac) <- "ATAC"
# VariableFeatures(atac) <- names(which(Matrix::rowSums(atac) > 50))
# atac <- RunLSI(atac, n = 30, scale.max = NULL)
# atac <- RunUMAP(atac, reduction = "lsi", dims = 1:30)
atac <- FindTopFeatures(atac, min.cutoff = 'q0')
# Seurat 3
atac <- RunLSI(atac, n = 30) 

## Warning: RunLSI is being moved to Signac. Equivalent functionality can be
## achieved via the Signac::RunTFIDF and Signac::RunSVD functions; for more
## information on Signac, please see https://github.com/timoast/Signac

## Warning: RunLSI is being moved to Signac. Equivalent functionality can be
## achieved via the Signac::RunTFIDF and Signac::RunSVD functions; for more
## information on Signac, please see https://github.com/timoast/Signac

## Warning: RunLSI is being moved to Signac. Equivalent functionality can be
## achieved via the Signac::RunTFIDF and Signac::RunSVD functions; for more
## information on Signac, please see https://github.com/timoast/Signac

## Performing TF-IDF normalization
## Running SVD on TF-IDF matrix
## Scaling cell embeddings

# Seurat 4
#atac <- RunTFIDF(atac)
#atac <- RunSVD(atac, n = 30, reduction.key = "lsi_", reduction.name = "lsi")
# drop first dim as often correlated with seq depth
atac <- RunUMAP(object = atac, reduction = 'lsi', dims = 2:30)

## Warning: The default method for RunUMAP has changed from calling Python UMAP via reticulate to the R-native UWOT using the cosine metric
## To use Python UMAP via reticulate, set umap.method to 'umap-learn' and metric to 'correlation'
## This message will be shown once per session

## 15:24:14 UMAP embedding parameters a = 0.9922 b = 1.112
## 15:24:14 Read 4076 rows and found 29 numeric columns
## 15:24:14 Using Annoy for neighbor search, n_neighbors = 30
## 15:24:14 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 15:24:15 Writing NN index file to temp file /var/folders/15/gs2v829x4gxcdv0tlgk_5l4m0000gn/T//RtmpFgCFIA/file4960660f015e
## 15:24:15 Searching Annoy index using 1 thread, search_k = 3000
## 15:24:16 Annoy recall = 100%
## 15:24:17 Commencing smooth kNN distance calibration using 1 thread
## 15:24:19 Initializing from normalized Laplacian + noise
## 15:24:19 Commencing optimization for 500 epochs, with 152690 positive edges
## 15:24:25 Optimization finished

# scRNA-seq import
rna <- CreateSeuratObject(counts=countHBC)
rna$tech <- '10x'
rna$celltype <- grpHBC
rna <- FindVariableFeatures(rna, selection.method = "vst", nfeatures = 2000)
rna <- ScaleData(rna)

## Centering and scaling data matrix

rna <- RunPCA(rna, npcs=20)

## PC_ 1 
## Positive:  Slc6a6, Rps27, Junb, Aqp3, Btg2, Fos, Fxyd3, Apoe, Krt5, Jun 
##     Gclm, Socs3, Ier2, Gpm6a, Ubc, Ctsl, Cyr61, Ifitm3, Laptm4a, Atf3 
##     Gstm2, Csrp1, Ptn, Pltp, Fosb, Dusp1, Igfbp2, Klf4, Krt14, Gas6 
## Negative:  Tmsb10, Cystm1, Atf5, Plaur, Adam8, Krt6a, Map1b, Tuba1a, Tubb2b, Cst6 
##     Ptma, Ecm1, Stmn1, Emp3, Uchl1, Upk3bl, Ceacam1, Clcf1, Tm4sf1, Gprc5a 
##     Cldn7, Tubb6, Lypd3, Cldn3, Lamc2, Phlda2, Ncam1, Cldn4, Mal, Sult2b1 
## PC_ 2 
## Positive:  Igfbp2, Cald1, Tmsb10, Ogn, Tmsb4x, Igfbp4, Adam8, Krt5, Apoe, Plaur 
##     Fam101b, Tspan6, Abca4, Sparc, Tm4sf1, Col23a1, Lamc2, Krt6a, Clcf1, Tmem108 
##     Slc29a1, Krt14, Tubb6, Ngf, Cst6, Oaf, Tubb2b, Ecm1, Hmga2, Car12 
## Negative:  Prdx6, Adh7, Gsta3, Sftpd, Cyp2f2, Ly6a, Tst, Slc16a11, Ifitm1, Reg3g 
##     Muc1, Por, Cbr2, Adh1, Lrrc26, Selenbp1, Mgst1, Gpx2, Sult1d1, Gstp1 
##     Cbr3, Cd14, Cyp4b1, Cxcl17, Gsto1, Serpinb11, Elf3, Ly6c1, Serpina3n, Epas1 
## PC_ 3 
## Positive:  Adh7, Aldh3a1, Selenbp1, Adh1, Ly6a, Oat, Defb1, Fmo2, Cbr3, Gstp1 
##     Tmem176a, Atp1b1, Tmem176b, Efemp1, Igfbp4, Upk1b, Sparc, Ly6c1, Tgm2, Igfbp2 
##     Mgp, Kcnj16, Sult1d1, Gmpr, Cdo1, Krt15, Nrtn, Cyp4b1, Muc1, Bsg 
## Negative:  Cldn4, X1810011O10Rik, Klf4, Sfn, Zfp36, S100a10, Mast4, Atf3, Tagln2, Fosb 
##     Avpi1, Arc, Nr4a1, Cldn7, Nfat5, Egr1, Klf6, Cdkn1a, Phlda1, Lmna 
##     Crym, Neat1, Anxa1, Krt18, S100a6, Acsm4, Lypd3, Papss2, Ier2, Krt23 
## PC_ 4 
## Positive:  mt.Co1, mt.Co3, mt.Atp6, Hpgd, Pla2r1, X1500009L16Rik, Acsm4, Rps27rt, Car12, Dapl1 
##     Myo6, Fth1, Papss2, Lypd2, Insl6, Cyb5r3, Scd2, Gas6, Trp63, Nmb 
##     Eya1, Piezo2, Limd2, Crym, Sema5a, Serpinb1b, Ftl1, X4631405K08Rik, Fndc1, Scgb1c1 
## Negative:  Actg1, Anxa3, Ier3, Nfkbia, F3, Tmem176b, Krt7, Sfn, Cyr61, Cldn4 
##     Krt14, Tacstd2, Tuba1c, Tmem176a, Fos, Bhlhe40, Junb, Tnfrsf12a, Gmpr, Cish 
##     Fst, Cxcl16, Plaur, Krt17, Ubc, Trib1, Oat, Atf3, Adam8, Dusp5 
## PC_ 5 
## Positive:  Fosb, Egr1, Neat1, Irs2, Pdcd4, Nr4a1, Lmna, Klf9, Cbx3, Dst 
##     Ccdc3, Mast4, Sdc4, Timp3, Nfia, Kdm6b, Birc5, Maff, Spry2, Il6st 
##     Stmn1, Hist1h2ap, Nfat5, Sik1, Arid5a, Gm26532, Pbk, Gclc, Tpm1, Fam101b 
## Negative:  Aldh2, Gstm1, Prdx1, Aqp3, Ces1d, Tst, Gsn, Clu, Ephx1, Cpe 
##     Ifitm3, Ubc, Krt18, Prdx6, Hspb1, Galm, Fmo6, Hspa5, Ppa1, Krt5 
##     Krt8, Pltp, Mgst1, Ctsl, Cyp2a5, Alas1, Cryl1, Calml4, Gstm2, Actg1

rna <- RunUMAP(rna, dims = 1:20, min_dist=.1)

## Warning: The following arguments are not used: min_dist

## 15:24:29 UMAP embedding parameters a = 0.9922 b = 1.112
## 15:24:29 Read 3064 rows and found 20 numeric columns
## 15:24:29 Using Annoy for neighbor search, n_neighbors = 30
## 15:24:29 Building Annoy index with metric = cosine, n_trees = 50
## 0%   10   20   30   40   50   60   70   80   90   100%
## [----|----|----|----|----|----|----|----|----|----|
## **************************************************|
## 15:24:29 Writing NN index file to temp file /var/folders/15/gs2v829x4gxcdv0tlgk_5l4m0000gn/T//RtmpFgCFIA/file49606322f6aa
## 15:24:29 Searching Annoy index using 1 thread, search_k = 3000
## 15:24:30 Annoy recall = 100%
## 15:24:31 Commencing smooth kNN distance calibration using 1 thread
## 15:24:32 Initializing from normalized Laplacian + noise
## 15:24:32 Commencing optimization for 500 epochs, with 122146 positive edges
## 15:24:37 Optimization finished

p1 <- DimPlot(atac, reduction = "umap", group.by="treatment", label=TRUE) + 
  NoLegend() + 
  ggtitle("scATAC-seq")

## Warning: Using `as.character()` on a quosure is deprecated as of rlang 0.3.0.
## Please use `as_label()` or `as_name()` instead.
## This warning is displayed once per session.

p2 <- DimPlot(rna, group.by = "celltype", label = TRUE, repel = TRUE) + NoLegend() + ggtitle("scRNA-seq")
p1 + p2

# Markers for activated/resting HBC in scRNA-seq
FeaturePlot(rna, features = "Krtdap")

FeaturePlot(rna, features = "Krt6a")

FeaturePlot(rna, features = "Krt5")

FeaturePlot(rna, features = "Krt14")

FeaturePlot(rna, features = "Trp63")

# Markers for activated/resting HBC in scATAC-seq
FeaturePlot(atac, features = "Krtdap") + xlim(c(-5,7))

## Warning: Could not find Krtdap in the default search locations, found in
## ACTIVITY assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 14 rows containing missing values (geom_point).

FeaturePlot(atac, features = "Krt6a") + xlim(c(-5,7))

## Warning: Could not find Krt6a in the default search locations, found in ACTIVITY
## assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 14 rows containing missing values (geom_point).

FeaturePlot(atac, features = "Krt5") + xlim(c(-6,6))

## Warning: Could not find Krt5 in the default search locations, found in ACTIVITY
## assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

FeaturePlot(atac, features = "Krt14") + xlim(c(-6,6))

## Warning: Could not find Krt14 in the default search locations, found in ACTIVITY
## assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

FeaturePlot(atac, features = "Trp63") + xlim(c(-6,6))

## Warning: Could not find Trp63 in the default search locations, found in ACTIVITY
## assay instead

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

## Warning: Removed 55 rows containing missing values (geom_point).

Assign sub-hybrid labels to RNA-seq

ctHybridSub <- as.character(rna$celltype)
hlpid <- which(rna@reductions$umap@cell.embeddings[,1] > 2 & rna@reductions$umap@cell.embeddings[,2]>-1)
ctHybridSub[hlpid] <- NA
hlp2id <- which(rna@reductions$umap@cell.embeddings[,1] > 2 & rna@reductions$umap@cell.embeddings[,2]>-1 & rna@reductions$umap@cell.embeddings[,2]<4)
ctHybridSub[hlp2id] <- "hybrid1"
hlp1id <- which(rna@reductions$umap@cell.embeddings[,1] > 2 & rna@reductions$umap@cell.embeddings[,2]>4)
ctHybridSub[hlp1id] <- "hybrid2"
rna$ctHybridSub <- ctHybridSub


DimPlot(rna, group.by = "ctHybridSub", label = FALSE, repel = TRUE) +
  ggtitle("scRNA-seq") 

table(ctHybridSub)

## ctHybridSub
## activated   hybrid1   hybrid2   resting 
##       602       130       148      2184

Explore hybrid1 markers from ATAC

expressionBoxplot <- function(rna, gene){
  counts <- as.vector(rna@assays$RNA[gene,])
  df <- data.frame(counts=log1p(counts),
                   celltype=rna$ctHybridSub)
  ggplot(df, aes(x=celltype, y=counts)) +
    geom_boxplot() +
    xlab("Cell state") +
    ylab("Log(Expression + 1)") +
    theme_bw() +
    theme(axis.text.x = element_text(angle=45, vjust=0.5)) +
    ggtitle(gene)
}

hybrid1Markers <- c("Muc5b", "Sult1c1")

FeaturePlot(rna, features = "Muc5b")

expressionBoxplot(rna, "Muc5b")

FeaturePlot(rna, features = "Sult1c1")

expressionBoxplot(rna, "Sult1c1")

Explore hybrid2 markers from ATAC

hybrid2Markers <- c("Pax3", "Pax7", "Pax9", "Dmrt2", "2610016A17Rik")

FeaturePlot(rna, features = "Pax3")

expressionBoxplot(rna, "Pax3")

FeaturePlot(rna, features = "Pax7")

expressionBoxplot(rna, "Pax7")

FeaturePlot(rna, features = "Pax9")

expressionBoxplot(rna, "Pax9")

FeaturePlot(rna, features = "Dmrt2")

expressionBoxplot(rna, "Dmrt2")